Formation of Polymorphic Cluster Phases for Purely Repulsive Soft Spheres 
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We present results from density functional theory and computer simulations that unambiguously 
predict the occurrence of first-order freezing transitions for a large class of ultrasoft model systems 
into cluster crystals. The clusters consist of fully overlapping particles and arise without the existence 
of attractive forces. The number of particles participating in a cluster scales linearly with density, 
therefore the crystals feature density-independent lattice constants. Clustering is accompanied by 
polymorphic bcc-fcc transitions, with fee being the stable phase at high densities. 

PACS numbers: 64.70.Dv, 82.30.Nr, 61.20.Ja, 82.70.Dd 



The distinguishing feature of soft matter systems is 
the vast separation of length and time scales character- 
izing the extent and motion of their constituent enti- 
ties. Whereas soft matter mixtures are typically solu- 
tions in a microscopic solvent, the solute particles are 
complex macromolecular aggregates of mesoscopic spa- 
tial dimensions ilj. The ability to control the architec- 
ture and chemical nature of these macromolecules, com- 
bined with the flexibility in influencing the solvent prop- 
erties and the composition of the system, gives rise to 
an unprecedented freedom in tuning the effective inter- 
actions between the particles and opens up the possi- 
bility to steer the macroscopic properties of the system 
P, ,2j. The richness of spontaneously forming complexes 
in soft matter encompasses length scales that exceed the 
dimensions of the individual macromolecules. Indeed, 
the latter can self-organize in a variety of ways, giving 
rise to so-called hypermolecular structures that en- 
compass a large number of mesoscopically-sized entities. 
Characteristic examples are the complex phases encoun- 
tered in ternary mixtures of oil, water, and amphiphilic 
surfactants or in block copolymer blends, as well as the 
emergence of cluster formation between colloidal parti- 
cles, which has attracted a grea t deal of attention re- 
cently H S B II II S EUl El- The underlying 
physical mechanism that drives the emergence of hyper- 
molecular structures is widely believed to rest on the exis- 
tence of competing interactions among the mesoscopic so- 
lute constituents. For example, the dominant mechanism 
that guarantees the stability of finite clusters in colloidal 
|M QilSl or biological |^ systems stems from the presence 
of short-range attractions and long-range repulsions in 
their effective interaction potential. Whereas the former 
provide the driving force for unlimited cluster growth, 
the latter act as a barrier against it. The efficiency of 
the barrier grows fast with increasing cluster population, 
therefore further accumulation of colloids into the clus- 
ters is brought to an end when a specific, optimal cluster 
occupancy is reached Cluster formation is a highly 

topical issue in current soft matter research, due to the 



large variety of cluster morphologies that form 0, ^ 
and to the relevance of these structures in influencing 
vitrification and gelation |^ |^ 0| . 

In this Letter, we report on a different mechanism that 
gives rise to a distinct type of cluster formation, and 
which does not rest on the explicit existence of compet- 
ing interactions. Contrary to the cases in Refs. 5,,^J^^, 
the constituent particles we consider are allowed to over- 
lap and are purely repulsive. Both conditions are readily 
fulfilled for various types of polymeric macromolecules, 
e.g., pol yme r chains [21 j polyelectrolytes 0|, or den- 
drimers [l5j |. Under certain, general conditions on the 
properties of the Fourier transform ^{q) of the interpar- 
ticle potential ^{r), we demonstrate that the particles 
form aggregates that further organize into regular clus- 
ter crystals with multiple site occupancy. We explicitly 
confirm the theoretical results by performing extensive 
computer simulations on a specific system that represents 
the entire class of effective interactions . 

Ultrasoft effective interactions 0, ^3 '^i'^^ many sur- 
prises in the topology of their phase diagrams and 
the types of crystal phases that arise, even for purely 
isotropic pair potentials <I>(r) Il9l l20l|. For the case 
in which the $(r) is non- negative and bounded, a gen- 
eral criterion determining the topology (but not the crys- 
tal structures) of the phase diagram has been put for- 
ward [T^ . If the Fourier transform $(<?) is non-negative 
(termed Q^-class), then the system displays re-entrant 
melting with an upper freezing temperature. If, on the 
other hand, ^(g) oscillates around zero (termed Q^- 
class), then a transition into an ordered cluster phase 
will take place at arbitrary temperatures. This corre- 
sponds exactly to the ordered 'clump phase' described 
in |23|. The argument put forward in Ref. rests 
on the fact that, except at small densities and tempera- 
tures, the fiuid state of the systems at hand is extremely 
well described by the mean-field approximation (MFA) 
c(r) = — /3$(r), where c(r) is the direct correlation func- 
tion and P = (fceT)"^, with Boltzmann's constant A;b 
and the absolute temperature T. Consequently, the fiuid 
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FIG. 1: The phase diagram of the GEM-4 model, as obtained 
by DPT. The shaded area represents the coexistence region of 
the liquid and the cluster bcc phase. The dash-dotted curve 
is the A-Iine of the system, calculated in the MFA. 



structure factor S{q) is given by S{q) = [1 + /3p^{q)]~^, 
where p denotes the number density. Consider now the 
Q^-class and let g, be the wavevector for which at- 
tains its minimum, negative value. Then, S{q) develops 
a real pole at along the so-called A-line 16, 22, 23|, 
which signals an instability of the uniform phase. Evi- 
dently, the A-line has the shape ksTx = i-e., it 
persists at all temperatures. Since g* is set solely by the 
interaction, the crystal lattice constant a (x q~^ should be 
density- independent, at least on the freezing line. Thus, 
the number of particles in the elementary cell should 
change accordingly, a requirement that can be fulfilled 
by the formation of multi-particle aggregates (clusters). 

Although the re-entrant melting scenario has been con- 
firmed [20L 0) the clustering scenario has re- 
ceived considerably less attention up to now. In this 
work, we explicitly demonstrate its validity. We per- 
form a detailed investigation of the generalized exponen- 
tial model of index n (GEM-n), <i>(r) = e exp[— (r/cr)"], 
with n — A. It can be shown that $(r) belongs to the Q"*"- 
class for n < 2 and to the Q^-class for n > 2. For n = 2, 
the Gaussian core model (GCM) of Stillinger [23| is re- 
covered. The GEM-4 is treated here as a representative 
of the Q^-class. Suitably tailored dendrimers that have 
been assembled in a computer simulation show evidence 
for a GEM- n- type of effective interaction with n > 2, 
hence this model reflects the behavior of realistic systems 

We define p* = pa^ and T* = k^T/e. Our investiga- 
tions consist of a combination of sophisticated minimiza- 
tion algorithms, density functional theory (DPT), and 
advanced Monte Carlo (MC) simulations. We started 
with a calculation at T = 0, allowing for the formation 
of clusters in which ric particles sit on top of each other, 
and minimized the lattice sum with respect to the crys- 
tal structure and ric. Since the periodic lattice that the 
clusters form is a priori unknown, we took advantage 
of an unbiased search strategy based on genetic algo- 
rithms [sj- Only fee- and bcc-crystals were predicted, 
which were used as candidates at finite temperatures. 
For this purpose, we employed DFT with the highly ac- 



curate MFA-excess free energy [23, 123, |29| J^^y\p] = 
(1/2) / dri / dr2/9(ri)p(r2)<i>(|ri - r2|). For the one par- 
ticle density, p(r), we made the Gaussian ansatz: p(r) = 
Pci(r-R.) = n,(a/7r)3/2 ^^^^^^ e-"('--R.)^ , where 
Pci(r) is the cluster density profile and the vectors 
{Ri} form a Bravais lattice. The total free energy is 
J^[p\ = -?id[p] + -^cx[p], with the ideal part ^id[p] = 
k^T J drp{r)[\n{p(r)A^) - 1], A being the thermal de 
Broglie wavelength. For aa^ > 20, Fjd can be approx- 
imated analytically and the variational free energy per 
particle, F/N = f{nc,a), takes the form: 



f{nc,cx) = k^T Innc -t- 3 ln(-\/c«T^77r) 4- 



-a{r+Rif/2 



$(r) 



+ {nc - 1)V fl / drr2e-"'-'/2$(r), 



(1) 



where the primed sum is carried over all lattice vectors 
excluding R = 0, Ri = |Ri|, and the irrelevant term 
3fcBrin(A/(T) has been dropped. The function f(ric,a) 
is then minimized at any state point with respect to ric 
and a. For a given crystal (fee or bcc) and density p* , the 
cluster population ric and the lattice vector lengths Ri are 
coupled to one another, Ri/a oc li{nc/ p*Y^^ , where h are 
lattice-specific numerical coefficients. At the minimum 
of f{ric,Oi), which corresponds to a mechanically stable 
crystal, the particular property Uc oc p* is fulfilled, so 
that the lattice constants of both the bcc- and the fcc- 
lattices remain fixed at all p*-values. In particular, the 
nearest neighbor distances for bcc and fee have the values 
rfbcc = 1.368(7 and dfcc = 1.414 cr respectively. 

In order to check also the thermodynamic stability of 
the crystals, we calculated the free energy of the uniform 
fluid employing the MFA, c(r) = —(3^{r), as a closure 
to the Ornstein-Zernike relation and following the en- 
ergy route to thermodynamics. With the free energies of 
the uniform and crystalline phases readily available, the 
phase diagram can be drawn, which is shown in Fig. ^ 
At low densities, the system is fluid. Upon compression, 
a first-order clustering transition occurs, leading at suf- 
ficiently high densities to the fcc-structure. Above the 
triple temperature, = 0.4, a wedge-shaped cluster 
bcc-region intervenes between the fluid and the cluster 
fee. Thus, the system also shows polymorphic transi- 
tions between cluster solids. The bcc-fcc density gap is 
very narrow, contrary to the fluid-crystal gap. In agree- 
ment with the arguments in Ref. If?!, the freezing and 
melting curves are almost straight lines and they precede 
the occurrence of the A-line in the phase diagram. 

To put the theoretical predictions in a stringent test, 
we also carried out MC simulations in the NVT-ensemble. 
Typically we used systems with approximately 5000 par- 
ticles and extended the computer experiments to 150 000 
MC-sweeps. Considerable speedup of the simulations was 
achieved by using a discretized simulation technique |32j | . 
The symmetry of the resulting regular structures was an- 



FIG. 2: Three simulation snapshots of a GEM-4 system for 
T* = 0.4 and p* = 2.5, 3.5, and 7 (from left to right). Particle 
diameters are not drawn to scale but are chosen to optimize 
the visibility of the structures. 




FIG. 3: Modified free energy density pPp* /N = pPp* /N + 
Kip* of the GEM-4 system against p* at T* = 1, as obtained 
by DFT and MC-simulations. A thermodynamically irrele- 
vant term K\p* has been added for clarity of presentation. 
Inset: a close-up on the region of stability of the cluster bcc- 
lattice. Here, a different linear term K2P* has been added. 
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FIG. 4: Cluster size, ric, and localization parameter, acr^, for 
the GEM-4 system at T* = 1, plotted against the density 
p*. Lines: DFT-results; points: MC-simulations. There are 
small discontinuities at the density of the bcc-fcc transition. 
The inset shows DFT-results that corroborate the ric oc p*- 
relation. The dashed part of the line is an extrapolation to 
low densities, for which the crystal is unstable. 
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FIG. 5: Cluster density profile, Pci{r), of the GEM-4 system 
at T* — 0.1 and p* — 9. Inset: Semi- logarithmic plot of Pci(r) 
against r^, demonstrating the Gaussian shape of the former. 



alyzed as proposed in Ref. '33i|. In the simulations, we 
found evidence of spontaneous clustering and crystalliza- 
tion, which is demonstrated by the snapshots in Fig. |21 
all for T* = 0.4. Whereas in the left panel (p* = 2.5) the 
system is obviously in the fluid phase, formed by isolated 
particles as well as clusters, at a higher density (p* = 3.5, 
middle panel), aggregates of particles have formed upon 
compression. These are still 'loose' and are located on 
a slightly distorted fee lattice. At a still higher density 
(p* = 7.0, right panel), the particles are tightly bound in 
clusters whose crystalline arrangement is evident. 

To determine the chemical potential p we used 
Widom's particle insertion [s^l supplemented by the 
overlapping distribution method |30ll34|. and the free en- 
ergies -Fbcc and Ffcc were obtained via the Gibbs-Duhem 
relation. A compilation of DFT- and MC-results for the 
free energies of all phases at T* = 1 is shown in Fig. 
13 It can be seen that the DFT-results are in excellent 
agreement with simulations, a fact that amply confirms 
the accuracy of the former and of the phase diagram in 
Fig. n In the inset of Fig. 13 we show a close-up of the 
three-phase region, to demonstrate that the cluster bcc- 
crystal is not preempted by a transition between fluid 
and cluster fee. 

The quantitative accuracy of DFT is not limited to 
thermodynamic quantities; it also holds for more de- 
tailed, structural ones. In particular, we have measured 



in the simulations the mean occupancy Uc of the clusters 
in the solid phases and compared it with ric obtained 
from DFT. Representative results for T* = 1 are shown 
in Fig.^ where the excellent agreement between the two 
can be seen. In particular, the linear dependence of Uc 
on p* is fully conflrmed. There is only a tiny difference 
in the population of the bcc- and fcc-clusters, see Fig. 0] 
This leads to to a ratio dfcc/^bcc — 1.034, very close to 
the value \^/\/3 = 1.029 obtained when ric is identi- 
cal for both lattices. At fixed p*, Uc is also practically 
T*-independent, so that a given crystal features a single 
lattice constant in the whole range of its stability. 

In addition, the clusters become more compact with 
increasing density: this is refiected in a corresponding 
increase of the a-values in the DFT approach. Indeed, 
as can be seen in Fig. the Gaussian density profile 
that was taken as an ansatz in the DFT approach is 
in excellent agreement with the one measured in MC. 
This allows for an extraction of a localization param- 
eter a from the simulation results, which can be com- 
pared to the one from DFT. The comparison is excel- 
lent at low temperatures and worsens somewhat as T* 
increases, however the discrepancies between the two re- 
main below 10% at T* = 1.0, as can be seen in Fig. 0] 
It appears that a also has a linear dependence on p*, a 
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feature unknown for crystals formed by harshly-repulsive 
particles. This property is intricately related to cluster 
formation and can be understood as follows. Let us fix 
all particles but one on their lattice sites and consider 
the potential V{r) they exert on this test particle, which 
performs small oscillations around a site at r = 0. Taking 
only the m nearest neighbor sites into account we have 
V{r) = {ric - l)$(r) + nc J2 - Rnn|), where the sum 
runs over all nearest neighbor clusters, located at the vec- 
tors Rnn of length d. For small r, V{r) can be expanded 
as V{r) ^ V{0) + {l/2)mnc<^"{dy, since $"(0) = [H. 
Thus, the test particle is a thermal oscillator and the den- 
sity profile p(r) cx exp[—(}mnc^" {d)r^ /2] results. The 
Gaussian shape is seen to arise from the local harmonic 
character of the site potential and the localization pa- 
rameter a can be identified as a = (3mnc^" (d)/2. Since 
Uc oc p* , the linear dependence of a on p* follows. 

All salient properties of the GEM-4-model that drive 
cluster formation are common to the entire Q^-class, 
thus the phenomena presented here should be general. 
The spontaneous formation of clusters appears counter- 
intuitive at first sight as, indeed, it occurs at the complete 
absence of competing interactions at the level of the pair 
potential. The underlying reason is rather the emergence 
of competing trends in the free energy, as can be seen in 
Eq. The entropy loss due to particle aggregation and 
the 'self-interaction' within the cluster (the kBTlnric- 
term and the last term on the rhs of Eq. respectively) 
disfavor the growth of Uc. However, increasing Uc also en- 
tails the gain in avoiding close contacts with the nearest 
neighbors, due to the concomitant increase of the lengths 
Ri in the second term on the rhs of Eq. At the same 
time, excessive growth of is also unfavorable, because 



it drives the inter-neighbor interaction term to zero, as 
the i?i's then lie outside the range of <i>(r): the entropic 
and self-interaction terms overtake as Uc — > oo and stop 
further aggregation. Interactions with the neighbors are 
also indispensable for the mechanical stability of the crys- 
tals, providing the required restoring forces for particle 
oscillations around the lattice sites. The stability of the 
clusters against both decomposition to Uc = 1 and un- 
limited growth towards ric — > oo is therefore seen to arise 
from a competition between m^ra-cluster interaction and 
entropy, on the one hand, and mier-cluster interaction, 
on the other. The necessary requirements for this sce- 
nario are that is bounded (to allow full overlaps) 
and that it falls rapidly enough to zero for r ~* oo, so 
that $(q) develops oscillations. 

We have presented a detailed analysis of the phase be- 
havior of a particular soft-interaction-system representa- 
tive of a broad class of effective interaction potentials 
that are realistic for ultrasoft, polymeric colloids. A 
novel mechanism for the development of cluster phases 
has been demonstrated to be at work, which gives rise 
to polymorphic crystals with unusual structural proper- 
ties. The system is accurately described by a mean- field 
density functional, as confirmed by comparison with com- 
puter simulations. Work along the lines of tailoring den- 
drimers that show Q^-interactions is under way. Future 
directions include the investigation of anomalous diffu- 
sion, gelation, and vitrification of such systems. 
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